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1.  Introduction 

In  the  course  of  supernova  remnant  (SNR)  evolution,  there  appear 

2 

to  be  four  well-defined  stages:  the  first  ~10  years  of  unimpeded 

expansion  of  hot  gas  away  from  the  collapsed  central  object;  a  period 
4 

~2xl0  years  during  which  this  envelope  sweeps  up  interstellar  material, 
slowing  down  somewhat  as  it  does  so;  a  period  ^XlO^  years  during 
which  line  emission  is  the  most  effective  cooling  mechanism;  and  a 
terminal  stage  in  which  the  expanding  material  becomes  indistinguishable 
from  the  background  gas. 

During  the  second  stage,  a  shock  wave  is  launched  which  propagates 
ahead  of  the  expanding  shell.  In  the  present  paper,  we  are  concerned  with 
the  hydrodynamic  character  of  this  shock  wave.  Other  workers  (see,  for 
example,  Chevalier  1976)  have  carried  out  elaborate  one-dimensional 
calculations  incorporating  realistic  transport,  radiation  and  atomic 
physics  models.  While  these  are  indispensable  for  acquiring  a  complete 
understanding  of  SNR  behavior,  idealized  simple  models  are  also  very 
useful.  They  afford  accurate  predictions  over  a  restricted  parameter 
range  and  qualitative  descriptions  over  a  wider  range. 

Such  a  model  has  been  derived  by  Sedov  (1946;  1959)  for  a  blast 
wave  expanding  from  a  point  source  into  a  gas-filled  surrounding  region. 
[An  abbreviated  description  has  been  given  by  Landau  and  Lifshitz  (1959); 
see  also  Newman  (1977).  Taylor  (1950) and  Von  Neumann  (1947)  independently 
developed  similar  models  in  studies  of  atmospheric  explosions.]  The  model 
consists  of  a  similarity  solution  of  the  ideal  gas  equations  for  a  strong 
shock  (Mach  number  M  =  V/cs  >>  1,  where  V  is  the  shock  speed  and  c 

s 

the  velocity  of  sound  in  the  undisturbed  medium)  produced  by  deposition  of 
Note:  Manuscript  submitted  December  20,  1979. 
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an  explosion  energy  W  at  the  origin  in  a  medium  of  density  p.  This 
model  can  be  expected  to  be  fairly  accurate  in  the  second  stage  of 
SNR  evolution,  during  which  the  shock  is  effectively  strong  and  line 
radiation  processes  negligible.  Chevalier  (1976)  and  others  Jo  in  fact 
utilize  it  in  this  way  in  studying  the  evolution  of  Type  II  supernovas. 

A  question  which  comes  readily  to  mind  is  whether  the  Sedov 

solutions  are  stable,  particularly  when  p  is  allowed  to  vary  with 

position.  In  their  general  form  the  solutions  are  quite  messy.  The 

fluid  variables  p,  v  and  p  are  given  implicitly  as  functions  of  r  and  t 

2-5  - 

through  the  similarity  variable  Wt  / p  r  (for  uniform  p ) .  Hence 
investigation  of  the  linear  stability  of  perturbations  about  these 
solutions  is  likely  to  be  impossible  in  the  general  case  except  through 
numerical  approximations . 

It  is  not  a  trivial  exercise  to  demonstrate  stability  even  of  plane 
shocks,  a  problem  which  has  been  treated  by  D'yakov  (1954)  andErpenbeck 
(1962)  in  ideal  gas-dynamic  systems,  and  by  Gardner  and  Kruskal (1964)  for 
MHD  shocks.  All  of  these  calculations  found  that  shock  waves  propagating 
in  a  uniform  medium  are  stable.  The  physical  mechanism  is  easy  to 
describe.  A  small  ripple  in  the  shock  front  gives  rise  to  divergence 
(convergence) in  the  curved  portion  ahead  (behind)  the  main  front.  These  regions 
become  weaker  (stronger)  than  the  unperturbed  shock,  hence  propagate 
slower  (faster)  than  average,  thus  reducing  the  amplitude  of  the  ripple. 
Evidently  the  longer  the  perturbation  wavelength,  the  weaker  the  stabilizing 
effect  of  this  mechanism. 

When  this  reasoning  is  applied  to  spherical  shock  waves  it  becomes 
obscure  and  must  be  regarded  as  no  better  than  a  plausibility  argument.  If 


2 


the  density  p  drops  sufficiently  rapidly  with  increasing  radius,  it 
is  possible  that  the  ripples  ahead  of  the  shock  run  away,  while  those 
behind  fall  farther  behind,  leading  to  instability.  Lerche  and  Vasyliunas 
(1976)  and  Isenberg  (1977)  have  claimed  that  for  at  least  some  decreasing 
power-law  density  distributions,  Sedov  shocks  are  unstable  for  any 
value  of  the  adiabatic  index  y.  Their  conclusions,  obtained  after  very 
elaborate  analysis,  were  remarkable  in  predicting  instability  at  short 
wavelengths,  contradicting  the  intuitive  argument  appropriate  to  the  planar 
case.  Newman  (1979)  has  disputed  the  results,  on  the  ground  that  Lerche  and 
Vasyliunas  (1976)  and  Isenberg  (1977)  improperly  treated  the  boundary 
condition  at  the  shock  front. 

The  present  paper  is  addressed  to  the  problem  of  determining  the 
stability  of  a  restricted  class  of  the  Sedov  solutions,  those  in  which 
the  outward  flow  behind  the  shock  in  homologous.  This  type  of  solution  was 
apparently  first  considered  by  Primakof f  (see  Courant  and  Friedrichs  1948) , 
who  was  studying  underwater  explosions.  He  found  that  for  p  =  const,  the 
ideal  fluid  equations  have  an  explicit  similarity  solution  when  y  =  7, 
the  approximate  value  of  the  effective  adiabatic  index  of  water  at  high 
pressure.  Keller  (1956)  generalized  the  results  to  arbitrary  power-law 
undisturbed  density  profiles  and  determined  the  value  of  the  density  exponent 
q  corresponding  to  each  choice  of  y  >  1. 

Using  a  formalism  originally  developed  by  Bernstein  and  Book  (1978)  and 
Book  and  Bernstein  (1979)  to  study  the  Rayleigh-Taylor  instability  of 
homologous  expansions  and  contractions,  we  solve  the  linearized  equations 
of  motion  exactly.  Our  criterion  of  stability  is  that  the  amplitude  of 
the  shock  front  perturbations  divided  by  the  shock  radius  vanish  as  t  +  *. 
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We  are  able  to  show  that  the  Primakoff  blast  wave  and  its  three-dimensional 
generalizations  are  stable  against  all  modes,  while  the  corresponding 
two-dimensional  line  blast  wave  solutions  are  stable  against  all  flute¬ 
like  (independent  of  z)  perturbations  (the  case  k ^  0  is  not  amenable  to 
treatment) , 

The  plan  of  the  paper  is  as  follows.  In  Section  2  we  review  the  derivation 
of  the  generalized  Primakoff  model,  using  an  approach  similar  to  that  of 
Keller  (1956).  In  Section  3  we  derive  the  linearized  equations  describing 
the  evolution  of  a  small  perturbation  about  the  basic  state.  Invoking  the 
boundary  conditions  obtained  from  the  Rankine-Hugoniot  relations  across  the 
shock  reduces  the  calculations  to  solution  of  an  eigenvalue  problem.  The 
results  are  presented  in  Section  A  for  cylindrical  and  spherical  blast  waves. 

In  Section  5  we  discuss  briefly  the  implications  of  our  results. 


2.  Generalized  Primakoff  blast  wave  model 


The  density  p,  velocity^  and  pressure  p  for  an  ideal  polytrope 
with  adiabatic  index  y  satisfy  the  equations 


p  +  p  V  *x=  °i  (2.1) 

P  JL+  Vp  =  0;  (2.2) 

P  +  YP  V  •  v.=  0,  (2.3) 


£ 

where  a  dot  denotes  the  material  derivative  —  +  v  ■  V.  On  a  surface 

O  t 

^Swith  normal  n,  moving  with  velocity  (where  both  ji  and JSL  can 
depend  on  position),  the  fluid  variables  can  change  discontinuously 
according  to  the  Rankine-Hugoniot  ("jump")  conditions 


^pn  •  0 L  -  jO)  =  0; 

(2.4) 

*  0L  ~  iy  -  jl)  +  P2^>  =  0; 

(2.5) 

\Cp  (V  -  v)2  +  n,  *  =  0, 

(2.6) 

where  <  >  denotes  the  jump  in  the  quantity  enclosed. 

By  dotting  and 

crossing  (2.5)  with we  find  the  projections 

,  p[n  •  (V  -  v)]2  +  p)  =0 

(2.7) 

and 

( P  n,  *  (V  -  v)  IjlX  (J,-v)])  =0. 

(2.8) 

The  latter,  by  virtue  of  (2.4),  reduces  to 

<  n  xg  -  v)>  =  0. 

(2.9) 

We  will  use  bars  to  distinguish  quantities  ahead  of  from 

those  behind.  Let  us  assume  that^  =  0,  i.e.,  the  fluid  in  front 

of  is  stationary.  It  follows  that  ji  x  jj,  =  0,  so  the  velocity 

behind  lies  in  the  direction  of  j^.  Hence  only  the  magnitudes 
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v  =  a.  •  j£  and  V  =  ja  •  enter  into  the  jump  conditions,  and  we 
can  rewrite  (2.4)  -  (2.6)  as 


p(V  -  v)  =  p 
P  (V-v)2+p=p 

P  (V  -  v)3  +  ^  (V  -  v) 

In  the  limit  of  very  strong  shocks, 
of  (2.10)  -  (2.12)  simplify  to 


v; 

(2.10) 

v2  +  p  ; 

(2.11) 

=  p  V3  +  V. 

(2.12) 

effectively  p  -*  0  and 

the  solutions 

P.  =  1+L 
P  Y-l* 


(2.13) 


v 

V 


"  Y+l  * 

2  -  „2 


p  =__p  v  =  p  v  V. 

Behind  the  shock  front  jS  we  seek  solutions  of  (2.1) 
with  homologous  (uniform)  expansion. 


(2.14) 

(2.15) 

(2.3)  consistent 


R  =  r  f (t) . 


(2.16) 


Here  R  is  the  position  of  a  particular  element  of  fluid  at  t,  and  r  is 
the  position  that  element  would  occupy  at  some  time  t^  (when  f  =  1) ; 
by  assumption  the  function  f  is  the  same  for  all  fluid  elements.  There 
is  no  special  physical  significance  to  the  choice  of  the  Lagrangian 
variable  r.  We  could,  for  example,  label  each  element  instead  by  the 
position  it  occupied  prior  to  being  overtaken  by  the  shock,  but  use  of 
r  simplifies  the  analysis. 

For  symmetric  one-dimensional  motion,  Eqs.  (2.1)  -  (2.3)  become 
p  +  PR1-V  |r(R^1v)  =  0;  (2.17) 
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(2.18) 


p  v  +  =  0; 

oR 

P  +  Y  p  R1_V  (RV_1v)  =  0,  (2.19) 

where  v=  1,  2,  3  for  planar,  cylindrical  and  spherical  geometry, 
respectively.  Using  the  Lagrangian  description  in  terms  of  r  and  t, 
we  find  from  (2.17) 


P (r , t)  =  PQ(r)f 


-  v 


(2.20) 


while  from  (2.19), 


p(r,t)  =  pQ(r)  f 


-  v  y 


(2.21) 


Substitution  into  (2.18)  now  yields  an  equation  which  separates  into 
a  spatial  and  a  temporal  o.d.e.,  viz., 


2 

o>  rp 


(2.22) 


and 


"  f  v  (y-1)+1 


2 

oi  . 


(2.23) 


Equation  (2.23)  can  be  integrated  twice,  yielding 


1  2  2m2  f  v (1-y) 

=  v (Y-l) 

and 

f  =  U  +  +  1]  [  y.  ]  15  w(t-t  ) 

2  J  u  v (y-l)  o 

From  (2.16)  the  velocity  satisfies 


(2.24) 


2 

v(y-l)+2  _  (2.25) 


v  =  R  =  rf  =  Rf/f.  (2.26) 

We  denote  the  radius  of  the  shock  at  time  t  by  S(t)  .  When  the 

jump  condition  (2.14)  is  applied  at  R  =  S,  we  obtain 

•  • 

Sf  _  2S_  , 
f  Y+l 
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(2.27) 


which  integrates  to 


1±I 


S  =  s  f 


where  s  is  the  shock  position  at  t  =  t^.  Hence  from  (2.13), 


p(S, t)  =  ^(S/f)f  =  p 


-v=  I±1  ' 

Y-l 


Suppose  p  has  a  pcwer-law  dependence  on  R.  In  that  case  we  can 


P  =  Pq  (R/s  )q  -  QR4, 


Pq,  Q  and  q  constant.  From  (2.28)  -  (2.30), 

.S,  Y+l  -  -  v  ,  S  .  q  Y+l  -  £  V+q  .  S  .  q 

Vf*  ■  ^  po£  (-r>  ■  ^  pof  <Tf> 


•  til;  <  _s_) 

Y-l  P0  V  s  f;  Y  , 


whence  for  arbitrary  r 


2+1  r  2 Ul+ai  + 

po(r)  =  p0  (~}  y-1 


Similarly,  from  (2.15), 


p(S,t)  -  P0(s/f)f"V1'.  ^  s2 


.  i!L  (  i  )1  s2  £?-i  f2  (r!i)2 
Y+l  's'  z  K  2  >• 


(2.28) 


(2.29) 
write 

(2.30) 


(2.31) 


(2.32) 


(2.33) 
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,  s,  =  ^0  2u^ _  Y+l  2  _S_  q  y-1+vy+v(1-y) 

PQ(  f  ’  Y+l  v(y-1)  (  2  )  (  s  '  * 

2  2 

(Y+l)  Prv  s  u  e  q  c  Y+  v-  1+q 

=  - y -  (-=-)  f 

v(y-1)  s  f' 


(Y+l)  Pq  s  2  w2 
v (y-1) 


q+  2 (y+ V-l+q) 

Hr7>  T'1 


(2  34) 


whence 


po(r)  = 


(Y+l)  Pq  S  2  W2 

V  (y-1) 


<^>,+ 


2 (y+ v -l+q) 


(2.35) 


Substitution  of  (2.32)  and  (2.35)  in  (2.22)  results  in  the  requirement 


Y(  v  -2)  -  3  v+2 
q  =  - - - 


(2.36) 


It  turns  out  to  be  convenient  to  write  all  equations  in  terms  of 


Y,  using  (2.36)  to  eliminate  q.  Thus  (2.32)  and  (2.35)  reduce  to 


,  \  Y+l  ~  /  r_\  v -2 
po(r)  =  7=1  po( 


tj  /  \  _  Y+l  -  2  2  ,  r  .  v 

Vr)  v(Y-l)  p0  s  “  (  s  )  » 


(2.37) 


(2.38) 


respectively.  Table  1  lists  the  values  of  q  for  two-  and  three- 
dimensional  flows  corresponding  to  some  typical  choices  of  y*  Note 
that  for  v=3  the  medium  ahead  of  the  shock  is  uniform  (q=0)  when 
Y=7,  the  original  Priraakoff  solution.  For  v=2,  however,  the 
undisturbed  medium  is  always  nonuniform  except  in  the  incompressible 
limit  y  +  “• 
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Table  1 


Y 

q 

v  =  2 

v=  3 

1 

-2 

-3 

5/3 

-3/2 

-2 

3 

-1 

-1 

7 

-1/2 

0 

oo 

0 

1 
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The  derivation  of  the  equations  describing  a  self-similar  blast 
wave  was  first  carried  out  in  terms  of  Lagrangian  variables  by  Keller 


=  2(y-1)P _ 

[v(y-D+2]2  t2 


by  (2.40)  and  (2.41). 

Equations  (2.41)  -  (2.43)  constitute  an  Eulerian  description 
of  the  flow  behind  the  shock.  Note  that  s  and  w  occur  only  in  the 

*(l-v) 

.  2-v  v(y-l)+2 

expression  s  w 

It  is  useful  to  rewrite  (2.41)  -  (2.43)  in  terms  of  Q  and  the 
total  energy  W, 

S 

W  «  ir2 V_  1  f  <*R  rV_1  (j  f^2  +  >  ,  (2‘44) 

*0 


instead  of  the  less  easily  interpreted  s  and  w.  For  this  purpose  we 
evaluate  (2.44)  [note  that,  by  (2.43),  the  two  terms  in  the  integrand 

are  equal],  with  the  result 


W  = 


Y+l  ~  2 

J j  Q  u 

(Y-l) 


2  [  v(y~1)+2  ] 
s  Y+l 


(2.45) 


Then  we  can  rewrite  (2.41)  as 


p(R.t) 


=  1+1 

"  Y-l 


v(Y-l)[v(Y-l)+2]2 
2  V  TT  (y+1) 


2  (1- v) 
v(y-1)+2 


v(y+1)  2(1-v) 

_  qv(y-1)+2  wv(Y-1)+2  R v- 2 


Ml-vJ _ 

t  v(y~1)+2  , 


(2.46) 


with  a  similar  expression  for  p  by  way  of  (2.43).  This  form  of  the 
solution  is  identical  with  that  given  by  Sedov  (1946;  1959). 
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3.  Linearized  equations  of  motion 

We  follow  Bernstein  and  Book  (1978)  in  obtaining  linearized 
equations  for  the  development  of  a  small  perturbation  about  the 
solutions  of  Section  2.  The  element  of  fluid  whose  trajectory 
is  given  byjl(r,t)  is  assumed  to  undergo  a  displacement  to 
R(r,t;  +  £>(r,t) .  The  perturbation  £  satisfies  the  linearized 

*0*  tv* 

form  of  (2.2) , 

4  +  Pii  +  -  V>  •  y  -  <3-» 

where  first-order  quantities  are  distinguished  by  the  subscript  1. 
Substituting  the  perturbel  velocity,  density  and  pressure  from 

vx  =  j,  ,  (3.2) 


and 

Px  =  -  ypvR  ’  jj_,  (3.4) 

and  making  use  of  the  relation  (2.16),  we  obtain 


-2  v(y-l)  +2  V  ,  _  ;  y  2-v 

a)  f  i+rV*b--Lr 

*"*’  —  *— ■  v 


V(r'V*£)  -  (V£)*r  =  0,  (3*5) 


where  V  denotes  the  gradient  with  respect  to  r. 

We  seek  solutions  of  (3.5)  by  the  method  of  separation  of 
variables.  Substituting 


i>(r,t)  =  X(r)T(t) , 


0.6) 
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we  find 


fv(y  l)+2  T  =  _  ^2t 


pX  -  rff  +  ^Lv(r2o)  +  (v-2)raj  +  V (r *X)  0, 


(3.7) 


(3.8) 


where  p  is  a  separation  constant  and  a  =  V*X.  Substituting  for 


f  from  (2.25)  converts  (3.7)  into 


+  —  >lT' 


(3.9) 


which  is  homogeneous  in  t.  Equation  (3.8)  is  likewise  homogeneous 
in  r  for  v  =  3,  and  also  for  v  =  2  if  X  is  independent  of  z,  the 
coordinate  in  the  direction  of  the  axis  of  symmetry.  It  follows 
that  r  and  t  both  enter  into  §>  with  power-law  dependences,  viz.. 


„  a 

X  ~  r  , 


T  ~  t^. 


(3.10) 


(3.11) 


By  (2.16),  the  solutions  expressed  in  terms  of  the  Eulerian  variables 
R, t  have  the  same  property  (with  different  coefficients  a',  6'). 


Writing 


X  =  ar  , 
r 


ro  =  br  , 


(3.12) 


(3.13) 


where  X^  is  the  radial  component  of  X,  we  see  that  (3.8)  yields 


(p+a)  a  +  L^(a+V“l)  -1 J  b  ■  0. 

v 


(3.14) 


A  second  equation  connecting  a  and  b  results  when  we  take  the 
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divergence  of  (3.8): 


£(a+l)  (a+2)  -  A]a  +  |y-l  +  [-^^-Y-lj  (vta-l) 

+  *  L(a+1) (a+2)  -  Aj|  b  =  0.  (3.15) 

Here  A  is  the  coefficient  of  the  terms  in  the  Laplacian  resulting  from 

the  angular  dependence  of  J>.  For  v  =  3,  A  =  £(£,+1),  while 
2 

for  v  =  2,  A  =  m  ,  where  l  and  m  have  their  usual  meanings. 

From  (3.9)  and  (3.11)  we  get  a  relation  between  y  and  6, 

=  _  Lv(Y-l)+2]2  B(e-l) 

P  2v(y-1)  *  (3.16) 

In  order  to  proceed  further,  we  need  to  apply  appropriate  boundary 
conditions  to  the  solution. 

As  a  result  of  the  perturbation,  the  shock  front undergoes 
a  small  displacement  into J  .  We  can  describe  this  by  saying  that 
at  time  t  a  point  on  iS  atJS^t)  is  shifted  to 

S' (t)  =  S  +  C(S,t) .  (3.17) 

«■»  «•> m 

There  is  an  arbitrariness  in  the  definition  of  L,  only  the  normal 
component  of  which  is  significant.  We  thus  free  to  define  0*  in  the 
most  convenient  manner,  as  a  mapping  along  the  unperturbed  normal  : 

C (S,t)  =  n  (S,t)  C(S,t)  (3.18) 

■»  ■» 

We  introduce  and  ,  the  first  order  corrections  to  ji  and  jy, 
respectively,  and  label  the  perturbed  fluid  variables  with  primes 
to  indicate  that  they  are  evaluated  on  A  instead  of  In  first  order, 
the  jump  conditions  (2.13)  -  (2.15)  yield 


0; 


(3.19) 


a -*i  +j>i  -i* < te. •  Ji  +.ci  •  s>;  <3-20> 

pi-m  te  -.v,(s  -ii+ii  •£>•  (3-21) 

From  kinematic  arguments,  n^  is  readily  shown  to  have  the  form 

n  =  n  n  •  (Vt)  •  n  -(VC)*  n,  (3.22) 

*»*  1  »•  * 


from  which  it  follows  that  n  •  =  0.  Thus  is  orthogonal  to  n, 

and  therefore  to  v  and  V.  Equations  (3,20)  and  (3.21)  simplify  to 


n  •  v,  = 

Ho.  «■! 


n  •  V. 

Ho  «*i 


4p  ^n  *  V) (n  •  V^) 

?  1  =  - ^+T -  =  2?  <“  *  V  &  '  z'l*  ' 


(3.23) 


(3.24) 


The  expressions  (3.2)  -  (3.4)  for  the  perturbed  fluid  variables 
must  now  be  replaced  by  their  Eulerian  counterparts. 


(3.25) 


Zi  -*“*•*  V5 


(3.26) 


Px  “  “  YP?R  *  ^  ~  J  *  VRP*  (3.27) 

Here  ^  (R, t)  is  the  displacement  in  position  experienced  by  a 
fluid  element  whose  unperturbed  position  at  time  t  was  R.  To 
evaluate  the  perturbed  quantities  at  S*,  we  expand  in  Taylor  series  to 
obtain 
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4  pSVV 

_  YPSZ  +  V p  U-Zr)  =  - —  +  qpc.  (3.36) 

By  eliminating  V  from  (3.34)  -  (3.36)  we  can  express  c  in  ferms  of  I , 

C  =  4 ( v-1)  SZ’  (3.37) 

and  then  substitute  for  ?  to  obtain  the  desired  boundary  condition 
at  R  =  S,  expressed  in  terms  of  jj  alone: 

(y-1)  si  +  [v(Y-l)  +  2]tcx  =  0.  (3.38) 

Since  by  (2.16)  RZ  =  re, we  can  use  (3.11)  -  (3.13)  to  write 
(3.38)  in  the  form 

(Y-l)b  +  ['Xy-D  +  2]8a  =  0  .  (3.39) 

Solving  (3.14),  (3.15)  and  (3.39)  for  a/b  yields 

(u+v- 1)  -  1 
a  v _ 

b  p+a 


h-1  + 


-  1]  (v+a-1)  [(a+1)  (a+2)-A ] 

(a+1) (a+2)  -  A 


Equations  (3.40),  together  with  (3.16),  constitute  a  set  of  three 
algebraic  relations  in  a,  S  and  p.  The  solutions  are  subject  to  an 
additional  condition,  namely  that  of  "regularity"  at  the  origin.  This 
is  the  requirement  that  the  second-order  contribution  to  the  total  energy. 


(3.41) 


/S 

dRRV_1{pi2  +  p[(y-l)(VR  •  Q‘ 


+  VrL  :  ?Ri]  } 


not  diverge  at  R  =  0.  Using  (2.16)  to  rewrite  the  integral  in  terms 
of  r,  we  see  that  it  is  convergent  provided 


Re  a  > 1  -  v. 


(3.42) 


which  is  also  the  condition  that  be  finite  at  r  =  0.  Thus  only  modes 
satisfying  (3.42)  are  physically  realizable. 


4.  Solution  of  the  eigenvalue  problem 

Equations  (3.16)  and  (3.40)  were  solved  for  both  v  =  2  and  v  =  3.  The 
degree  of  this  system  is  not  easy  to  determine  by  inspection,  and  initially 
solutions  were  obtained  numerically  for  various  choices  of  A.  When 
it  became  apparent  that  there  are  only  four  solutions,  the  equations  were 
manipulated  to  reduce  them  to  a  quartic  in  a,  g  or  y  (these  forms  and 
their  solutions  were  obtained  largely  by  use  of  the  interactive  computerized 
symbolic  manipulation  system  MACSYMA) . 

As  our  criterion  of  stability,  we  evaluated  the  time  dependence  of 
£/S.  If  Q/S  increases  with  time,  the  basic  state  is  unstable;  otherwise 
it  is  stable.  From  (2.28)  and  (3.11)  we  see  that  this  ratio  is 
proportional  to 


(S/f)a  t^  f  %[(Y-l)a-y-l]  B 

(Y+D/2  ~  1  C 

f 

Hence  the  condition  for  stability  is 

Ref  0, 


where 


r  = 


(y— 1)  ot-Y-1 
v(Y-l)  +  2 


(4.1) 


(4.2) 


(A. 3) 


Accordingly  we  list  below  the  values  of  i  corresponding  to  those 
found  for  u  and  6. 

(a)  v  =  2 

The  four  roots  of  the  a  equation  can  be  expressed  in  the  form 


(y* 


1)m  =  1  ± 


D©E‘ 


(4.4) 
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where 


D  =  [y2  -  (Y2-1)A  f2 


(4.5) 


and 

E  =  y2(Y2-1)(1-A)  +  2(2-y2)  ±  D.  (4.6) 

Here  the  circled  sign  in  (4.4)  varies  independently  of  the  uncircled 
ones  in  (4.4)  and  (4.6),  which  are  either  both  (+)  or  both  (-).  The 
corresponding  values  of  8  are  given  by 

2y(Y+De  =  Y(1  ±  D^E*5,  (4.7) 

while  from  (4.3)  those  of  I  are 


(4.8) 


The  associated  values  of  u,  which  we  omit,  can  also  be  obtained  as  solutions 
of  a  quartic  or  by  substitution  in  (3.16). 

When  y  4  1  or  when  A  =  0  or  A  =  1,  Dis  real.  The  latter  cases 


correspond  to  m  =  0  and  m  =  1,  the  two  flutelike  modes  with  the 


longest  wavelengths  allowed  in  cylindrical  geometry,  and  are  therefore 
the  ones  expected  (according  to  the  simple  picture  discussed  in  Section  1) 


to  be  most  unstable. 
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Table  2 


G- 


0 


0 


.  1±L_ 
2y 


Table  2  lists  the  associated  values  of  a,  8  and  i  for  each  of  the  four 
roots,  labeled  by  the  circled  and  uncircled  signs.  We  note  that  (3.42) 
is  satisfied  for  all  solutions  found  in  these  limits  except  the 
root  when  m  =  0.  Furthermore,  F  ^  0  for  all  cases  shown  in  the  Table. 

For  sufficiently  large  values  of  y  and  A,  D  becomes  imaginary  and 
the  real  part  of  T  always  equals  -1/2.  Thus  we  have  as  a  general  conclusion 
F  £  0,  so  that  when  v  =  2  the  blast  waves  solution  under  consideration 
are  stable  against  perturbations  with  any  value  of  m. 

(b)  v-  3 

The  analysis  of  the  spherical  case  exactly  parallels  that  carried 
out  for  v  =  2.  The  four  a  roots  are  given  by 


2 (y2-l)a  =  4-y-y2  ±  D^±)e5, 


(4.9) 


where  now 


D  =  [(3y-D2  -  4(y2-1)A]1s 


(4.10) 


E  =  yA(9-4A)  +  6y3  -  y2  (26-4A)  -  18y  +  37 


±  2(6-y-3y  )D, 


(4.11) 


while  the  corresponding  8  roots  are  given  by 


2(y+1)(3y-i)S  =  3Y-1  ±  yn(rph. 


(4.12) 


where  the  same  conventions  as  before  are  employed  with  regard  to  circled 


23 


and  uncircled  signs.  Substitution  of  (A. 9)  and  (4.12)  in  (4.3)  yields 


r  = 


D 

2 (3y-l) 


(4.13) 


When  y  -*■  1,  D  ->  2.  For  £  =  A  =  0,  D  =  3y  -  1,  while  for  £  =  1, 

A  =  2  and  D  =  3  -  y.  Table  3  shown  the  values  of  a,  B  and  T  in  these 
limiting  cases,  where  for  conciseness  we  have  written 


F  =  y4  +  12y3  -  34y2  -  36y  +  73.  (4.14) 

As  before,  (3.42)  holds  in  all  these  cases  except  the  ((^)-)  mode  for 
£  =  0,  and  T  £  0  for  all  cases  shown  in  the  Table.  From  (4.13),  it  is 
clear  that  the  latter  conclusion  holds  in  the  general  case  as  well, 
although  (3.42)  is  usually  satisfied  only  for  one  pair  of  roots. 


Table  3 


5.  Conclusion 


The  results  presented  in  this  paper  establish  rigorously 
the  stability  of  a  restricted  class  of  the  general  Sedov  solutions. 

This  class (the  Primakoff  blast  wave  and  its  analogs)  has  been  shown 
by  Sedov  (1959)  to  be  degenerate.  That  is,  the  solution  collapses  to 
a  single  point  in  the  phase  plane  defined  by  the  reduced  flow  and  sound 
speeds,  as  a  consequence  of  the  boundary  condition  being  applied  at  a 
singularity (in  this  case,  the  origin).  Oppenheim  et  al.  (1972)  make 
clear  the  relationship  of  such  degenerate  solutions  to  the  general  case. 

The  simplifying  assumption  responsible  for  the  tractability  of  the 
problem  we  have  solved  is  the  requirement  (2.36)  that  the  unshocked  medium 
have,  for  any  given  value  of  y,  a  particular  power-law  density  profile. 

As  far  as  we  know,  there  are  no  observations  of  supernovas  for  which 
the  envelope  gas  density  profiles  are  sufficiently  accurately  known  to  say 
whether  they  are  close  to  those  of  the  present  model.  It  would  be 
surprising  if  anything  so  simple  actually  occurred. 

The  results  presented  here  provide  no  basis  for  any  conclusions 
regarding  the  stability  of  the  class  of  Sedov  solutions  as  a  whole.  Never¬ 
theless,  it  is  probable  that  other  profiles  which  are  qualitatively 
similar,  i.e.,  decreasing  as  some  other  power  of  the  radius  or  perhaps 
logarithmically,  are  likewise  stable.  The  reason  for  saying  this  is 

that  (4.13)  predicts  positive (better  than  marginal)  stability  for  almost 
all  modes.  Even  a  hypothetical  shift  in  the  direction  of  instability  could 
well  leave  ReT  <  0  for  all  modes,  or  all  except  l  ■  0. 

We  note  that  the  analysis  of  the  cylindrical  problem  is  similarly 
incomplete.  It  is  incomplete  in  another  way  because  of  the  omission  of 
modes  with  k^  f  0.  Although  there  are  few  astrophysical  examples  of 
cylindrical  explosions,  they  do  arise  in  laboratory  experiments,  such  as 
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those  relating  to  the  propagation  of  intense  pulsed  charged-particle 
beams  through  gas-filled  chambers.  There,  too,  the  question  of  the 
stability  of  the  resulting  shock  waves  propagating  into  the  surrounding 
medium  has  been  widely  discussed.  In  such  experiments,  however,  the 
usual  case  of  interest  is  that  of  radially  increasing  density  profiles 
created  as  a  result  of  channeling  or  hole-boring  by  previous  pulses. 

Note  that  our  analysis  applies,  mutatis  mutandis,  to  imploding 
shocks  for  which  the  flow  behind  (i.e.,  outside)  the  front  is  homologous. 
This  is  a  restricted  case  of  the  general  solution  found  by  Guderley 
(1942).  If  we  replace  t  by  -t,  all  equations  remain  correct  except 
those  defining  the  energy,  where  the  integrals  must  now  be  carried  out 
from  S  to  «>.  Instead  of  (2.44)  we  find  a  divergent  expression  for  W. 
Equation  (3.41)  for  can  also  diverge  but  must  do  so  less  rapidly. 

The  condition  for  this,  which  is  also  the  condition  _hat  first-order 
quantities  be  small  compared  with  their  unperturbed  counterparts,  is 
readily  seen  to  be  a  i  1,  replacing  (3.42).  This  means  that  a  different 
subset  of  the  solutions  found  in  Section  4  corresponds  to  physically 
realizable  modes.  The  general  conclusion  ReT  £  0  implies  instability,  since 
C/S  diverges  as  -t  -*  0.  From  Table  3  we  see  that  the  fastest  growth  (T  =  -1) 
is  found  as  y  -*■  1  or  when  l  =  0,  but  we  conclude  from  (4.13)  that  insta¬ 
bility  persists  even  when  l  •*■  <*>. 

Considerably  more  work  can  probably  be  done  along  the  lines  of 
that  presented  here.  While  an  analysis  in  closed  form  of  the  stability 
of  the  general  Sedov  solutions  is  unlikely,  it  would  be  surprising  if 
the  use  of  more  powerful  mathematical  techniques  than  those  employed 
here  could  not  widen  the  range  of  accessible  cases. 
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